
#pragma once
#include "FEM2D.h"
/*!
 * \class CFEM2DVector
 *
 * \brief FEM class for vector unknowns.
 *
 * \author wanyzh
 * \date Jan 2021
 */
class CFEM2DVector :
	public CFEM2D
{
public:
	CFEM2DVector(int n=1);
	~CFEM2DVector();

public:
	void setNumVec(int n);
	int getNumVec();

	///@brief Initialize
	///
	///This function is inherited from the base class CFEM2D. It is used to initialize the memory of the boundary nodes and edges according to CFEM2DVector::m_nSizeVec
	void Init();

	///@brief Generate initial vectors
	Rsv generate_Initial(double(*fun)(CPoint),vector<CPoint>Pb);

	///@brief A general assemble function integration over 2D domain, see also CFEM2D::assembleMatrix2D.
	///The different between this version with the version in the base class CFEM2D::assembleMatrix2D is that 
	///the Pb,Tb, and basis function are the parameters. 
	///This is for the mixed FEM whose basis functions and finite element mesh are different for different unknowns.
	///@param basis_trial	Basis trial function.	
	///@Tb_trial			Tb for trial function.
	///@param basis_test	Basis test function.
	///@param Tb_test		Tb for test function.
	///@param matrix_size	The matrix size get from Pb, it is a pair and the first element represents the rows' number and second one is the cols' number.	
	///Other parameters and return are the same with the function CFEM2D::assembleMatrix2D.
	SpMat assembleMatrix2D(double(*coe)(CPoint, double t), double t,
		CBasisFunction &basis_trial, const int &basis_trial_der_x, const int &basis_trial_der_y, vector<CElement>&Tb_trial, 
		CBasisFunction &basis_test, const int &basis_test_der_x, const int &basis_test_der_y, vector<CElement>&Tb_test, const pair<int,int> &matrix_size);

	///@brief Assemble right side vector
	///
	///This function is simliar with the CFEM2D::assembleVector2D. The different is that the basis function and Tb are the parameters.
	///@param basis_test	Basis test function.
	///@param Tb_test		Tb for test function.
	///Other parameters and return are the same with the function CFEM2D::assembleVector2D.
	Rsv assembleVector2D(double(*load)(CPoint, double t), double t, 
		CBasisFunction & basis_test, const int &basis_der_x_test, const int &basis_der_y_test, vector<CElement> &Tb_test);
	

	///@brief Override the boundary conditions treatment function for vector unknown
	///it is related to the m_nSizeVec
	void treatBoundaryConditions(double(*coe)(CPoint, double t), double t, SpMat &A, Rsv &b,vector<unsigned int> num_of_values);

	///@brief stack matrix horizontally.
	///reference parameters to avoid the copy of the object.
	SpMat matrix_hstack(const SpMat &m1, const SpMat &m2);
	SpMat matrix_hstack(const SpMat &m1, const SpMat &m2, const SpMat &m3);

	///@brief Stack matrix vertically.
	SpMat matrix_vstack(const SpMat  &m1, const SpMat &m2);
	SpMat matrix_vstack(const SpMat &m1, const SpMat &m2,  const SpMat&m3);


	///@brief Stack vector
	Rsv vector_stack(const Rsv  &b1,const Rsv  &b2);
	Rsv vector_stack(const Rsv  &b1,const Rsv  &b2, Rsv const &b3);

	///@brief export result to vtk
	void export2VTK(string filename, const Rsv &x);
	void export2VTK(string filename, const Rsv &x, vector<CPoint> &P, vector<CElement> &T);

	///@brief Assign boundary edges value (r and value) for the index-th boundary which is the order defined in the .neu mesh file.
	///
	///This function call the CMesh::assignBoundEdgesValue to implenment.
	void assignBoundEdgesValue(unsigned int index, int vec_index, BOUNDARY_TYPE type, double(*value)(CPoint, double t), double(*r)(CPoint, double t));

	///@brief Assign boundary node value (r and value) for the index-th boundary which is the order defined in the .neu mesh file
	void assignBoundNodesValue(unsigned int index, int vec_index, BOUNDARY_TYPE type, double(*value)(CPoint, double t));

	///get mesh object CMesh.
	CMesh &getMesh();
	
	vector<CElement>& getTbTrial();
	vector<CElement>& getTbTest();
	vector<CPoint>&   getPbTrial();
	vector<CPoint>&  getPbTest();

	vector<CElement>& getT();
	vector<CPoint>&  getP();

	///get the matrix A's condition number 
	double getMatCondNum(SpMat A);

	///@brief impose a fixed value for the node.
	///
	///Sometimes we need to set the value of only one node and the node is not in the boundary_nodes, so we have to handle this  separately.
	///@param index		The index of the finite element node needed to set the value.
	///@param value		The value.
	void imposeFixValue(const int &index, const double &value,  SpMat &A,  Rsv &b);

	///@brief Solve linear equations  with MKL parallel method.
	///@param num_pro The number of the parallel processors.
	Rsv solveLinearEquaParallel(const SpMat  &A,const Rsv   &b, const int &num_pro);
	

protected:
	///@brief Gauss Quadrature over 2D element.
	///
	///This function is the similar with CFEM2D::GaussQuad2DCoeTrialTest.
	///@param basis_trial	The basis trial function.
	///@param basis_test	The basis test function.
	double GaussQuad2DCoeTrialTest(double(*coe)(CPoint, double t), double t,const vector<double> &Gauss_weights, const vector<CPoint> &Gauss_nodes, const vector<CPoint> &vertices,
		 CBasisFunction &basis_trial, const int &basis_index_trial, const int &basis_der_x_trial, const int &basis_der_y_trial,
		 CBasisFunction &basis_test, const int &basis_index_test, const int &basis_der_x_test, const int &basis_der_y_test);


	///@brief Gauss quadrature over 2D element.
	///
	///This function is similar with CFEM2D::GaussQuad2DLoadTrial.
	///@param basis_test	The basis test function.	
	double GaussQuad2DLoadTrial(double(*load)(CPoint, double t), double t, const vector<double> &Gauss_weights,const vector<CPoint> &Gauss_nodes,const vector<CPoint> &vertices,
		CBasisFunction &basis_test, const int &basis_index_test,const int &basis_der_x_test,const int &basis_der_y_test);

protected:
		
	int m_nSizeVec;///<Size of the unknown vector. m_nSizeVec=2 in 2D elasticity problem and =2 in mixed FEM for Stokes equations.
};

